
\chapter{关于 c++, 以及其他工具}

\section{c++ }
国内理科学生一般都要在本科阶段学习 c 语言，或者 c++。
我们这个课程基本都是做一些简单的面向过程的函数式编程，所以对我们来说，c 或者 c++ 都可以。
为了帮助大家温习，我们简单地回顾一下 c++ 的最基本内容。
如果有需要，可以在后面的学习过程中翻阅c++教科书。
国内使用比较广泛的是\cite{谭浩强c++}，国外比较受欢迎的有\cite{C++Primer}(也有中文翻译版，虽然我觉得原版更流畅)。

\subsection{安装 Dev C++}
我估计大家都用的 windows 操作系统，如果没有安装 c++ 编译器的话，我建议大家使用 dev-c++，因为它非常轻便简洁，也满足我们的需求。
安装和运行过程如下：
\begin{itemize}
	\item 下载 Dev C++ 安装文件，我附上下载地址\footnote{https://sourceforge.net/projects/dev-cpp/}。
	\item 下载完成以后，一般双击这个文件即可，但有时可能需要管理员权限，则右键这个文件，选择以管理员权限运行。
	\item 点击两三次 "OK" 或者 “下一步”即可完成安装。
	\item 建立新文件：File -> New -> Source File，然后编写程序，并保存为你指定的文件名。
	\item 编译：Execute -> Compile，或者快捷键 F9，界面下方会显示编译结果，包括错误的个数，警告的个数（警告的内容不一定是错误的，有可能只是编译器疑惑了）。
	\item 运行：Execute -> Run, 或者快捷键 F10，也可以用“编译\&运行”，快捷键为 F11。
\end{itemize}

\subsection{编译，运行：Hello world}
下面是最简单的 c++ 程序
\begin{lstlisting}
#include<iostream>   // 包括标准输入输出库
using namespace std; // 用命名空间
// 如果不写这一行，我们想用输入输出命令，就不能只打 cin, cout，而要写 std::cin, std::cout

int main(){
//主函数，包括函数返回类型（int)，主函数名（main），参数（无参数），和函数体（花括号内的内容）
	cout<<"Hello world!"<<endl;
//输出“Hello world!"这些字符。endl 是换行。
	return 0;
//主函数返回值为0，表示正常结束。你可以试试改成-1，然后看看运行结果。
}
\end{lstlisting}

在 Dev-C++ 中新建文件，编写好上面的程序以后，就可以用界面上的快捷键编译，运行。
如果窗口一闪而过，可以包括头文件stdlib.h，在 return 0;的前面加一行 system("pause")，让程序在 return 0 之前暂停。
system(char*)函数包括在cstdlib头文件中，如在windows下，就暂停进程，执行windows系统命令char*，如在linux下，就执行shell命令char*。
\begin{lstlisting}
#include<iostream>   
#include<cstdlib>
using namespace std; 

int main(){
    cout<<"Hello world!"<<endl;
    system("pause");
    return 0;
}
\end{lstlisting}
这样，你就可以看见屏幕上出现 Hello world! 的字样。
现在人们教任何编程语言，第一个都教这个程序，因为它仿佛一个初生的小机器人，在向全世界打招呼。

如果你输入的程序有错误，编译器会提示你哪里有错。
编译器反馈的信息不一定准确（毕竟它不是人类），但一般能反映出错误的大致位置。
然后你就按照提示，把输错的地方改正，直到编译通过。
如果你和李逵一样粗心，可能会需要多一点时间。。。
但是所有人类都是会犯错的，即使有多年编程经验的人，也自然会打错字。
It's human to make mistakes!

下面，我们就介绍得稍微快一点了，我们将回顾数据类型、数组、逻辑和循环语句、函数式编程、指针与动态内存，基本上就够我们用了。

\subsection{常用的数据类型}

所有语言的基础都是变量，只是有的需要特别声明变量的数据类型，有的不需要（比如python）。
在 c++ 中，大家经常用到的数据类型如下

\begin{table} [h]
	\centering
\begin{tabular}{|c|c|c|}
	\hline \hline
	数据类型 & 内存（byte) & 数值范围 \\
	\hline
	char	& 1			& $-128 \sim 127$ \\
	bool	& 1			& $true / false $ \\
	int		& 2 		& $-2^{31} \sim 2^{31}-1$ \\
	long int & 4		& $-2147483648 \sim 2147483647$ \\
	float	& 4		& $10^{-44} \sim 10^{38}$ \\
	double 	& 8		& $10^{-322} \sim 10^{308}$\\
	\hline \hline
\end{tabular} 
\caption[data-type]{本讲义常用的 c++ 数据类型。}
\end{table}
byte 是 “字节”，我们平时说的 kb, mb, gb, tb 中的 b 都是指 byte。
1 个 byte = 8 个 bit，bit 是“比特”，即二进制的一位。

c++ 中有加减乘除等等运算符，这些运算符有优先级规则，可以查书，也可以百度。
下面是一个简单的加法程序
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
	int a, b, c;// 声明三个整型变量
	cout<<" Please input values of integer a and b."<<endl;// 请输入两个整数
	cin>>a;//输入 a,b 的值
	cin>>b;
	c = a + b;// c 的值为 a 和 b 的和
	cout<<" The sum is : "<<c<<endl;
	return 0;
}
\end{lstlisting}
编译通过以后运行，屏幕上会出现提示：Please input values of integer a and b.
然后输入 3，回车，5，回车，它就会告诉你：" The sum is : 8"。
这就是一个很简单的整数加法计算器了。

\paragraph{数组}
如果我们需要同时操作 10 个整数，就可以定义一个 长度为 10 的数组，如下
\begin{lstlisting}
	int k[10];
\end{lstlisting}
也可以直接赋初始值 1,2，..., 10
\begin{lstlisting}
	int k[10]={1,2,3,4,5,6,7,8,9,10};
\end{lstlisting}
注意，在 c/c++ 中，数组的下标是从 0 开始的，如果要改变数组 k 中的第一个数字的值，使它变为 100，应该写
\begin{lstlisting}
	k[0]=100;
\end{lstlisting}

类似地，也可以定义浮点数数组，字符数组。

\subsection{逻辑语句 if, 循环语句 for 和 while}
\paragraph{逻辑语句 if}
下面我们写一个程序，输入两个双精度浮点数，让电脑告诉我们哪个更大
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
	double a,b;
	cout<<" Please input two float numbers:"<<endl;
	cout<<"a=";
	cin>>a;
	cout<<"b=";
	cin>>b;
	if(a>b)cout<<" Conclusion: a > b "<<endl;
	else cout<<" Conclusion: a <= b "<<endl;
	return 0;
}
\end{lstlisting}
其中就用到了 if 语句。
if 语句的格式是这样的
\begin{lstlisting}
    if(test){
        body1
	}
	else{
		body2
	}
\end{lstlisting}
其中 test 往往是一个逻辑语句，比如上面程序中的 a>b，那么 a>b 这个语句要么是真的，要么是假的，如果是真的，就执行 body1 中的所有语句，否则，就执行 body2 中的所有语句。

当然，还可以写得更复杂
\begin{lstlisting}
	if(test1){
		body1
	}
	else if(test2){
		body2
	}
	else{
		body3
	}
\end{lstlisting}
可以加上任意层的 else if，来区分各种不同的情况。

\paragraph{循环语句 for}
大家可能都听过小学生高斯的故事，关于“高斯求和”的技巧。
我们可以写一个程序，做 1 到 100 的总和。
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
	int i,sum=0;// 声明两个整型变量，其中 sum 的值为0
	for(i=1;i<=100;i++){
		sum += i;
	}
	cout<<" The sum of 1,2,...,100 is: "<<sum<<endl;
	// 输出结果
	return 0;
}
\end{lstlisting}
for 语句的结构是这样的：
\begin{lstlisting}
for(initial; test; increment){
	body;
}
\end{lstlisting}
其中，initial, test, increment, body 都是语句。
计算机看到 for 语句，将会做以下的事情
\begin{itemize}
	\item[1] 执行 initial，在上面的程序第 6 行中，就是令 i=1，即赋初值。
	\item[2] 判断 test 语句是真或假，即 i<=100 是否成立
	\item[3] 如果 test 语句是真的，就执行 body 中的所有语句,即第 7 行，让 sum 的值加上 i 的值。然后执行 increment 语句，即 i++，表示 i 的值加上 1，其值变为 2。然后回到第 2 步。
	\item[4] 如果 test 语句是假的，就结束这段 for 语句的工作。
\end{itemize}
所以，在上面的程序中，第 7 行会被执行 100 次，这就叫做循环。

\paragraph{死循环}
如果我们将上面的代码改成下面这样
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
	int i,sum=0;// 声明两个整型变量，其中 sum 的值为0
	for(i=1;1>0;i++){
		sum += i;
	}
	cout<<" The sum of 1,2,...,100 is: "<<sum<<endl;
	// 输出结果
	return 0;
}
\end{lstlisting}
即把 test 语句写成 1>0，那么 1>0 是永远成立的，所以程序将会无限次重复步骤 2 和 3，除非系统或者人为干预，程序永远不会停下来。
这就叫做死循环。


\paragraph{循环语句 while}
从 1 加到 100，也可以用 while 语句写：
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
	int i=1,sum=0;// 声明两个整型变量，其中 i 的值为1， sum 的值为0
	while(i<=100){
		sum += i;
		i++;
	}

	cout<<" The sum of 1,2,...,100 is: "<<sum<<endl;
	return 0;
}
\end{lstlisting}
while 的结构如下
\begin{lstlisting}
while(test){
	body
}
\end{lstlisting}
其中 test 和 body 都是语句，计算机看到 while 语句，会做如下事情：
\begin{itemize}
	\item [1] 检验 test 是否为真
	\item [2] 如果 test 为真，则执行 body 中的所有语句，然后回到步骤 1
	\item [3] 如果 test 为假，则结束整个 while 语句
\end{itemize}
你可以试试，把这个 while 语句也写成死循环。

\paragraph{变量的作用域}
c++中，需要注意变量的作用域。
上面的例子中，变量都是定义在函数体内部，那么，它们的作用域就是这块函数体，出了函数体以后，就无法调用这些变量了。
如果在一个for循环中定义变量，其作用域就是循环体，出了for循环，这个变量的内存就会被释放掉，无法再调用它。
所以下面的程序会报错。
\begin{lstlisting}
#include<iostream>
using namespace std;

int main(){
    for(int i=0;i<1;i++)
        i;
    cout<<"i="<<i<<endl;
    return 0;
}
\end{lstlisting}

\paragraph{算法：冒泡法排序}
所谓算法，就是一系列有序的动作，完成一个特定的任务。

下面是我的青椒炒肉丝攻略：先去超市买菜，再切青椒、肉丝，用酱油、料酒、姜丝、盐腌肉丝一刻钟，大火放油炒青椒（加少量水）到青椒稍微有一点软，大火放油炒肉丝，肉丝都发白时加入青椒，继续炒几分钟然后放生抽（上色）和少量盐，炒匀以后装盘放上餐桌。

这就是青椒炒肉丝的算法，它定义了一系列有序的动作（买菜、备料、炒菜等一系列细节），来完成特定的任务（做出青椒炒肉丝这盘菜）。

下面我们接触第一个数值算法，即冒泡法排序。
我们的任务是，给定任意 5 个不同的整数，写一个程序对它们进行自动排序，并按从小到大的顺序输出。
冒泡法的步骤如下
\begin{itemize}
	\item [1] 比较 $a_1, a_2$，如果 $a_1 > a_2$，就将它们调换。
	\item [2] 比较 $a_2, a_3$, 如果 $a_2 > a_3$，如果 $a_2 > a_3$ 就将他们调换。
	注意这里的 $a_2$ 只是指第 2 个数，它不一定等于步骤 1 中的 $a_2$ 了。
	\item [3] ...
	\item [4] 比较 $a_4, a_5$，如果不是从小到大，就调换。
	\item [5] 重复步骤 1-3
	\item [6] 重复步骤 1-2
	\item [7] 重复步骤 1
\end{itemize}
步骤 1-4 实现了一件事：将最大的数放到最末尾；
步骤 5 则将第二大的数放到次末尾；
步骤 6 将第三大的数放到了倒数第三位；
步骤 7 将第四大的数放到了倒数第四位；
任务完成，输出结果。

这些步骤，可以用两层嵌套的 for 语句来实现：
\begin{lstlisting}
#include<iostream>
using namespace std;
int main(){
	int a[5] = { 1,8,9,5,3};//待排序的数字
	int i,j,k;// 声明变量
	for(i=4;i>=1;i--){
		for(j=0;j<i;j++){ // 将 a[0] ~ a[i] 中的最大值，放到第 i+1 个位置 
			if(a[j]>a[j+1]){
				k=a[j]; // 交换 a[j] 和 a[j+1]
				a[j]=a[j+1];
				a[j+1]=k;
			}
		}
	}
	cout<<" The five numbers after sorting: "<<endl;
	for(i=0;i<5;i++){
		cout<<a[i]<<",";
	}
	cout<<endl;
	return 0;
}
\end{lstlisting}

\subsection{函数式编程}
如果工作比较复杂，我们在 main 函数中写几百上千行，就显得非常臃肿。
我们常常需要写一些功能函数，让它们互相调用，将 main 函数作为最外层的调用。
这种编程方式就叫做函数式编程，也是本课程使用的编程方式。
下面我们写一个比大小的函数，然后用 main 函数调用它。
\begin{lstlisting}
#include<iostream>
using namespace std;

void compare(int a, int b){ // 函数不返回任何值，所以类型为 void
	cout<<" Conclusion: ";
	if(a>b)cout<<a << ">" << b <<endl;
	else if(a==b)cout<<a << "=" << b <<endl;
	else cout<< a << "<" <<b <<endl;
}

int main(){
	int a,b;
	cout<<" Please input two integers:"<<endl;
	cout<<" a=";
	cin>>a;
	cout<<" b=";
	cin>>b;
	compare(a,b);//调用 compare 函数
	return 0;
}
\end{lstlisting}

如果定义的函数太多，我们可能难以找到 main 函数，所以，我们希望在 main 函数的后面定义 compare 函数。
c++ 要求所有函数被调用之前都被声明，所以我们需要在 main 函数前面加上 compare 的声明，然后在 main 函数的后面定义 compare 函数，如下
\begin{lstlisting}
#include<iostream>
using namespace std;

void compare(int a, int b);//声明 compare 函数

int main(){
	int a,b;
	cout<<" Please input two integers:"<<endl;
	cout<<" a=";
	cin>>a;
	cout<<" b=";
	cin>>b;
	compare(a,b);//调用 compare 函数
	return 0;
}

void compare(int a, int b){ // 函数不返回任何值，所以类型为 void
	cout<<" Conclusion: ";
	if(a>b)cout<<a << ">" << b <<endl;
	else if(a==b)cout<<a << "=" << b <<endl;
	else cout<< a << "<" <<b <<endl;
}
\end{lstlisting}

\subsection{递归函数}
一个函数可以调用自己，这样的函数叫做递归函数。
下面的代码是一个例子
\begin{lstlisting}
double factorial(int n){
    if(n<0){
        cout<<"error: invalid input for factorial().\n";
        //报错，\n是换行符，所以输出这一行以后屏幕上的光标会换行。
        exit(1);//异常退出，结束程序
    }
    else if(n==0) return 1;// 0!=1
    else return n*factorial(n-1);// n! = n* (n-1)!
}
\end{lstlisting}
在上面的代码中，我们定义了一个函数叫做“阶乘”，如果计算factorial(5)，它会返还5*factorial(4)，factorial(4)还是未知的，所以程序会继续求解factorial(4)，又返还4*factorial(3)，factorial(3)是未知的，于是程序自动求解下去，直到遇到factorial(0)=1，然后捡起前面的因子1*2*3*4*5，得到factorial(5)=5!。

递归函数的优点是，代码非常简洁，可读性强，代码往往是一个递推公式，在上面的代码中是$n! = n*(n-1)!$，通过这个递推公式让程序自动寻找解；缺点是递归深度是有限的，而且有时候内存的局域性不是很好。

原则上，任何递归函数都可以写成循环，但对于简单的任务，我们可能不在乎递归的上述缺点，而享受它代码的简洁清晰。
递归层数多了以后，函数头的开销会比较多，但c++的编译选项中有-O2和-O3，可以进行尾递归优化，但这是另一个话题了，等我有心情了，可能把尾递归优化的内容放到这里。。。

\subsection{指针与动态内存}
\paragraph{指针}
指针是一种特别的数据类型，它存放的是内存的物理地址。
请看下面的例子
\begin{lstlisting}
int a=3;
int * b;
b = & a;
cout<<" originally a="<<a<<endl;
*b +=3;
cout<<" later a="<<a<<endl;
\end{lstlisting}
上面的代码中，定义了整型变量a，其初值为3，整型指针（可以定义不同类型的指针）b。
\&是取地址符，所以第三行是用a的内存地址给b赋值，也就是说，让指针b指向变量a。
*b表示指针b指向的内存单元，所以第5行让那块内存单元的值增加3。
所以输出的结果为
\begin{lstlisting}
 originally a=3
 later a=6
\end{lstlisting}

\paragraph{动态内存}
内存的分配，分为两种；一种是编译阶段就决定了分配多大的内存，比如前面的所有程序；另一种是运行时才决定分配多大的内存。
后者叫做动态内存，用法如下
\begin{lstlisting}
int n;
cout<<"n="; cin>>n;//注意这里n必须是正整数，否则下一行将会报错
int *a = new int [n];
delete [] a;
\end{lstlisting}
第一行定义了正整数n，第二行输入n的值，第三行根据n的值，自动开辟n个整数的内存，后面可以进行计算，但最终别忘了像第4行一样，释放动态内存。
如果不释放动态内存，那些内存单元将无法重复利用，知道整个程序结束，我国人民将这样的内存叫做“僵尸内存”。

下面就是一个用动态内存的程序，自动给任意多个数排序。
\begin{lstlisting}
#include<iostream>
using namespace std;
int main(){
        int n;
        int i,j,k;
        cout<<" How many integers are to be sorted?\n n=";
        cin>>n;
        int *a = new int [n];//开辟动态内存
        for(i=0;i<n;i++){
                cout<<" a["<<i<<"]="; cin>>a[i];
        }
        for(i=n-1;i>=1;i--){
                for(j=0;j<i;j++){ // 将 a[0] ~ a[i] 中的最大值，放到第 i+1 个位>置 
                        if(a[j]>a[j+1]){
                                k=a[j]; // 交换 a[j] 和 a[j+1]
                                a[j]=a[j+1];
                                a[j+1]=k;
                        }
                }
        }
        cout<<" The "<<n<<" numbers after sorting: "<<endl;
        for(i=0;i<n;i++){
                cout<<a[i]<<",";
        }
        cout<<endl;
        delete [] a;//释放动态内存
        return 0;
}
\end{lstlisting}

如果你需要用动态内存定义多维数组，可以如下编写
\begin{lstlisting}
#include<iostream>
using namespace std;
int main(){
        int n;
        double **A = new double * [n];　//A为指向指针的指针
        for(int i=0;i<n;i++){
            A[i] = new double [n]; //A[i]为指针
            for(int j=0;j<n;j++){//给第i行赋值
                if(j!=i){
                    A[i][j]=1;
                }
                else
                    A[i][j]=0;
            }
        }
       
        for(int i=0;i<n;i++){//释放动态内存
            delete [] A[i];
        }
        delete [] A;
        
        return 0;
}
\end{lstlisting}
类似地，可以定义$n$维数组。

在上面定义的矩阵A中，其矩阵元的内存地址顺序为：
\begin{equation}
A[0][0], \cdots A[0][n-1], A[1][0], A[1][1], \cdots, A[1][n-1], \cdots
\end{equation}
利用这个信息，我们可以让程序少做内存的跨越，让程序更高效。

\section{用 gnuplot 自动作图}
gnuplot是一个非常简单的画图软件。
它可以用命令行进行交互式的调整（即打几行命令，就可以让它自动重画，直到你满意为止），所以熟练以后操作很快。
它也可以非常方便地与latex结合起来，结合gnuplot的画图优势与latex对公式编辑排版的优势。
它还可以制作3d图、动画gif文件，我们将做简单演示。

\subsection{gnuplot安装运行，第一张gnuplot图片}
\paragraph{安装运行}
在 ubuntu 系统下，输入
\begin{lstlisting}
sudo apt-get install gnuplot
\end{lstlisting}
然后选择 “Y”，即可安装 gnuplot。
在命令行中输入
\begin{lstlisting}
gnuplot
\end{lstlisting}
即出现 gnuplot 的提示符
\begin{lstlisting}
gnuplot>
\end{lstlisting}

在　windows　系统下，可以免费下载gnuplot安装文件\footnote{http://gnuplot.sourceforge.net/download.html}。
然后一路按提示安装即可。
安装完成以后，桌面会出现一个gnuplot图标，叫做“gnuplot..."，双击就会弹出gnuplot窗口，并出现gnuplot提示符
\begin{lstlisting}
gnuplot>
\end{lstlisting}

\paragraph{第一张图片}
我们可以在提示符下输入如下内容并回车，
\begin{lstlisting}
gnuplot>plot x*x with line
\end{lstlisting}
就会弹出一张图片，画着抛物线$y=x^2$。
plot和splot是gnuplot的基本命令，后面的x*x是函数，也可以改成sin(x),cos(x),exp(x)等等，gnuplot内置了这些常用的函数，幂次可以用$x**2,x**3$等表达。
with line是指用曲线作图，你可以试试with point。

我们可以让它显示x、y坐标轴标记，并修改曲线标记：
\begin{lstlisting}
gnuplot>set xlabel "x"
gnuplot>set ylabel "y"
gnuplot>plot x*x with line title "f(x)=x^2"
\end{lstlisting}
set xlabel "xxx"即将x轴名字标为xxx，set ylabel类似，而title后面双引号内则是曲线的标记内容。

在窗口的上端，有"export plot to file"的快捷键，存为png文件即可。
我们也可以改变输出方式，让它自动生成一个文件，如下
\begin{lstlisting}
gnuplot>set terminal eps
gnuplot>set output "1stGraph.eps"
gnuplot>set xlabel "x"
gnuplot>set ylabel "y"
gnuplot>plot x*x with line title "f(x)=x^2"
\end{lstlisting}
terminal即终端，在这里表示输出的方式，terminal默认的是wxt，即弹出图片。
现在我们把它设为eps，即生成eps文件，但我们还需要set output命令来指定文件名。
如果在Windows下操作，gnuplot的默认工作目录是“我的文档”，所以数据文件与生成的图片都在“我的文档”中。
可以通过界面左上角的“change directory...”改变工作目录。

我们也可以将这些命令写成一个脚本，并自动生成png、eps、ps等文件。
新建文本文档"1stGraph.txt"，其中内容为
\begin{lstlisting}
set terminal eps
set output "1stGraph.eps"
set xlabel "x"
set ylabel "y"
plot x*x with line title "f(x)=x^2"
\end{lstlisting}
然后在gnuplot窗口中打入
\begin{lstlisting}
gnuplot>load "1stGraph.txt"
\end{lstlisting}
就会自动执行1stGraph.txt脚本，生成1stGraph.eps文件。

我们可以在latex中包括这个eps文件，latex代码为
\begin{lstlisting}{language=tex}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\linewidth]{1stGraph}
\caption{gnuplot:第一张图}
\label{fig:1stGraph}
\end{figure}
\end{lstlisting}
效果如图\ref{fig:1stGraph}。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\linewidth]{1stGraph}
\caption{gnuplot:第一张图}
\label{fig:1stGraph}
\end{figure}

\subsection{散点图、折线图}
如果有一个呈两列的数据文件"ScatterDiagram\_LineChart.txt"，其中内容如下
\begin{lstlisting}
-10     2
-8      3
-6      1
-4      -1
-2      4
\end{lstlisting}
那么我们可以写如下脚本
\begin{lstlisting}
set terminal cairolatex size 4 inch, 2.4 inch
set output "ScatterDiagram_LineChart.tex"
set xlabel 'log$_{10}(h)$'
set ylabel '$\epsilon = \int^\infty_1 f(x) dx $'
plot "ScatterDiagram_LineChart.txt" title "折线图" w l, "ScatterDiagram_LineChart.txt" title "散点图" w p　pt 7
\end{lstlisting}
其中cairolatex是terminal的一种选项，它会生成tex文件，其中包含了图片与文字部分，文字部分用latex编辑，所以可以自由地编辑公式。
size 4inch, 2.4inch　表示生成图片的尺寸为4英寸x2.4英寸，这可以自由调节。
在set xlabel，set ylabel中我们展示了latex公式编码。
plot命令中，w l是with line的缩写。
plot命令后面可以跟两条画图命令，用逗号隔开，可以画一张散点图、一张折线图并显示在同一张图中。
pt 7　是散点图的一个选项，表示pointtype 7，即用内置的第7种点作图，如果想知道内置的点都是什么样的，可以在gnuplot中打test，就会弹出一张图片，给出相关信息。

注意，如果gnuplot的运行目录，与数据文件不在同一个地方，可能会报错。
这时可以使用gnuplot的上方第四个小图标“change directory”，找到数据文件所在文件夹，即可。

用gnuplot运行上述脚本以后，可以在latex中包括tex文件，如下生成图片
\begin{lstlisting}{language=tex}
\begin{figure}[h]
\centering
\input{ScatterDiagram_LineChart.tex}
\caption{gnuplot:散点图与折线图，x轴与y轴的图例没有任何意义，只是为了展示latex的公式编辑。}
\label{fig:ScatterDiagram_LineChart}
\end{figure}
\end{lstlisting}
效果如图(\ref{fig:ScatterDiagram_LineChart})。
\begin{figure}[h]
\centering
%\includegraphics[width=0.7\linewidth]{ScatterDiagram_LineChart}
\input{ScatterDiagram_LineChart.tex}
\caption{gnuplot:散点图与折线图，x轴与y轴的图例没有任何意义，只是为了展示latex的公式编辑。}
\label{fig:ScatterDiagram_LineChart}
\end{figure}

\paragraph{在图片中添加文字、箭头、直线}
\begin{lstlisting}
set label "hahaha" at -8,1
set arrow from -7.5,1.5 to -8,2.5 lw 5 linecolor rgbcolor "red"
set arrow from -8,0.8 to -6.6,0.8 nohead lw 2 linecolor rgbcolor "blue"
\end{lstlisting}
上面的三行分别在图上添加一段文字、一个红色箭头和一段蓝色横线，效果如图(\ref{fig:ScatterDiagram_LineChart2})所示。
\begin{figure}[h]
\centering
\input{ScatterDiagram_LineChart2.tex}
\caption{gnuplot:在图片中增加文字、箭头、直线。}
\label{fig:ScatterDiagram_LineChart2}
\end{figure}

\paragraph{在png等term的文字中添加希腊字母}
如果使用了latex相关的terminal，添加希腊字母是非常容易的，但是如果使用png等terminal，可以如下添加希腊字母
\begin{lstlisting}
{/Symbol abcdefghikl}
\end{lstlisting}
则会显示$\alpha\beta\chi\delta\epsilon\phi\gamma\eta\iota\kappa\lambda$，如果需要其他希腊字母，可以试出来，也可以随手百度。

\subsection{gif文件：动图}
gnuplot的terminal可以有很多选项，其中包括gif，并包括animate即动画选项。
下面我们介绍一个最简单的动画：波形图。
下面是gnuplot脚本
\begin{lstlisting}
set terminal gif animate
set output "sinx.gif"
set xrange [0:4*pi]
set yrange [-1:1]
do for [i=0:80]{
        plot sin(x+i*pi/40) linewidth 3 title sprintf("t=%d s",i)
}
set output
\end{lstlisting}
其中set xrange约定x轴范围，而do for循环生成81桢，这81桢自动生成gif文件。

用latex生成的pdf展示gif文件有点麻烦，需要将每一帧保存为图片，然后用animate包中的animategraphics命令生成，并且只有一些支持javascript的pdf阅读器才能显示，所以我就不在这个讲义中展示了。
后面的一些计算例子会生成动图，

\subsection{多图}
gnuplot可以将多张图在一起，下面是一个非常简单的示例脚本
\begin{lstlisting}
set multiplot layout 2,2    #　4张图以 2x2　排列
plot x                      #第一张
plot x**2                   #第二张
plot x**3                   #第三张
plot x**4                   #第四张
unset multiplot             #结束多图模式
\end{lstlisting}
在gnuplot提示符后面粘贴上面的代码，按下回车键，就会弹出图\ref{fig:multiplot}，
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{multiplot}
\caption{gnuplot:多图演示}
\label{fig:multiplot}
\end{figure}

我们还可以将两张图重叠放起来，比如一张大图里有一张小图
\begin{lstlisting}
set multiplot
plot x*x            #大图
set origin 0.3, 0.5 #设置小图原点位置（相对于大图）
set size 0.4, 0.4   #设置小图大小（相对于大图）
replot              #重复 plot x*x，即小图和大图是同一根曲线
unset multiplot     #结束多图模式
\end{lstlisting}
效果如图\ref{fig:multiplot2}
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{multiplot2}
\caption{gnuplot:多图演示}
\label{fig:multiplot2}
\end{figure}

\subsection{参数图}
如果想画圆，可以用参数坐标
\begin{eqnarray}
x &=& \cos t,\\
y &=& \sin t.
\end{eqnarray}
以t为变量，自动画出一个圆。
用如下脚本
\begin{lstlisting}
set terminal cairolatex size 7 cm, 7 cm
set output "circle.tex"
set parametric
set dummy t
set trange [-pi:pi]
plot cos(t), sin(t) with line linewidth 3
\end{lstlisting}
即可得到如下图形
\begin{figure}[h]
\centering
\input{circle.tex}
\caption{gnuplot:参数图。}
\label{figure:circle}
\end{figure}

\subsection{三维图}
plot命令用于画二维图，而splot可以画三维图。
比如，我们可以画出$f(x,y)=10x^2 + y^3$的曲面图，gnuplot脚本如下
\begin{lstlisting}
set terminal cairolatex size 5 inch, 4 inch
set output "3d_surface.tex"
set xlabel "x"
set ylabel "y"
set zlabel "z"
f(x,y)=10*x**2 + y**3
set isosamples 100
set view 60,60
splot f(x,y) title "$z=10x^2+y^3$"
\end{lstlisting}
其中，f(x,y)=10*x**2 + y**3定义了一个函数，也就是说gnuplot可以自由定义函数。
set view 60,60是指，观看角度如此定下：绕x轴旋转60度，z轴旋转60度。
效果如图\ref{figure:3d_surface}。
\begin{figure}[h]
\centering
\input{3d_surface.tex}
\caption{三维曲面图：$f(x,y)=10x^2 + y^3$}
\label{figure:3d_surface}
\end{figure}

上面是用解析函数做出曲面图，我们也可以用离散的数据点做出三维图。
作为示例，我们用一段c++代码生成一个曲面上的坐标点，然后用pm3d画出这个曲面。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#define grid 200

int main(){

        int i,j;
        double x,y,z;
        ofstream fp("test.dat");
        for(i=1;i<=grid;i++){
                x = -1 + (i-1)/(double)grid * 2;            // x的范围是[-1,1]
                for(j=1;j<=grid;j++){
                        y = -1 + (j-1)/(double)grid * 2;    // y的范围是[-1,1]
                        z = x*x + y*y;                      // z = 
                        fp<<x<<"\t"<<y<<"\t"<<z<<"\n";
                }
                fp<<"\n"; // 这一行很重要，x相同的数据块之间，空出一行
        }
        fp.close();
        return 0;
}
\end{lstlisting}
注意：x相同的数据块之间，空出一行。
运行上面的程序以后，得到数据文件test.dat，然后用gnuplot作图
\begin{lstlisting}
gnuplot> splot "test.dat" with pm3d
\end{lstlisting}
得到图片
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{data_to_3d}
\caption{离散的数据点做三维图。}
\label{figure:data_to_3d}
\end{figure}

\subsection{上色、灯光}
如果做的更加复杂一点，可以画出三维几何体，并给它上色，打上灯光。
下面的例子来自官网demo。
\begin{lstlisting}
set terminal cairolatex color size 4, 3
set output '3d_lighting.tex'
unset border
unset key       #unset the title
#set style fill   solid 1.00 noborder # set style fill density = 1.00
set dummy u, v  #independent variables
set object  1 rect from screen 0, 0 to screen 1, 1
set object  1 behind lw 1.0  dashtype solid fc  rgb "gray"  fillstyle   solid 1.00 border lt -1 #set background color-filling
set parametric # from normal functions to parametric functions
set view 236, 339, 1.245, 1 #set the view angle: rot_x, rot_z, scale, scale_z
set isosamples 75, 75
unset xtics
unset ytics
unset ztics
set title "PM3D surfaces with specular highlighting"
set urange [-pi,pi]
set vrange [-pi,pi]
set pm3d depthorder
set pm3d lighting primary 0.5 specular 0.6
set palette rgbformulae 8, 9, 7
splot cos(u)+.5*cos(u)*cos(v),sin(u)+.5*sin(u)*cos(v),.5*sin(v) with pm3d,     1+cos(u)+.5*cos(u)*cos(v),.5*sin(v),sin(u)+.5*sin(u)*cos(v) with pm3d
\end{lstlisting}
unset border, unset xtics等等，表示取消边框、xyz轴数字，unset key表示取消图例。\\
set dummy u,v表示u,v为自有参数。\\
set object两行定义的是正方形灰色背景。\\
set isosamples 75,75决定的是线的稠密度，如果稠密度比较低，画出来的更像多面体。\\
set pm3d depthorder表示离观察者更近的线条会挡住后面的线条，这样才会产生立体感。\\
set pm3d lighting primary 0.5 specular 0.6表示打上灯光，光强为0.5（最大值为1），而高亮效果的强度为0.6（最大值为0.6）。\\
set palette rgbformulae 8,9,7定制调色板，即确定颜色填充方案，如果改动数字，会得到不同的颜色填充效果，但是不是任何方案都好看。。。\\
最后一行splot将两个参数图画在同一张图上，即图\ref{figure:3d_lighting}中的两个环状体。
\begin{figure}[h]
\centering
\input{3d_lighting.tex}
\caption{三维图:背景为灰色色块，两个环状体由两个参数图实现，上色、灯光由palette, pm3d lighing实现。}
\label{figure:3d_lighting}
\end{figure}

在gnuplot界面中打入help，即可调出help文档，可以根据关键字或者目录查阅说明，探索更细致的功能。
本讲义中用的最多的是折线图、散点图，所以就介绍到这里了。

\section{openmp: 共享内存式并行}
\subsection{openmp: hello world}
openmp 是一个非常简单易上手的并行工具，在下文中我们将展示，只需要加两三行，就可以实现并行。
在 openmp 中，不同的线程可以共享一些内存（即公有一部分变量或数组等），所以它是一种共享内存式并行工具。
openmp 可以靠共享的一些变量实现不同进程的协作。

在后文中，我们将简单介绍 OpenMPI，它是一种分布式内存并行工具，即不同的进程不共享任何数据。
不同的进程向彼此发送信息来实现协作。

在 openmp 的官网上，有一份比较清晰的介绍课件
\footnote{https://www.openmp.org/wp-content/uploads/Intro\_To\_OpenMP\_Mattson.pdf},
以及相关的练习答案
\footnote{https://www.openmp.org/wp-content/uploads/Mattson\_OMP\_exercises.zip},
如果有兴趣，可以前往下载阅读。

现在的 c++ 编译器中，一般都已经带了 openmp 的库文件，只需要调用。
在 Dev-c++ 中，需做如下设置：工具栏->工具->编译选项->勾选“编译时加入如下命令”->输入“-fopenmp”。
然后我们介绍 openmp 的 hello world 程序
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<omp.h>//包括 openmp 的库函数，其中有对 omp_get_thread_num() 等等函数的声明。
			//c++ 要求每个函数在被引用之前都被声明

int main(){

	#pragma omp parallel// 创造新线程，每个线程重复执行花括号中的内容
	{
		int ID = omp_get_thread_num();//获取线程编号（第一个线程的编号是 0）
		cout<<" Thread "<<ID<<": hello world"<<endl;
	}

	return 0;
}
\end{lstlisting}
运行以后，会得到好多个 "Thread n: hello world"，如果你的电脑有 8 个线程，就会出现 8 个 hello world, n 为 0 到 7。
另外，你可能注意到，这些输出是乱糟糟的，这是因为这些线程独立地运行自己的任务，可能一个说出 hello，还没来得及说出 world，另一个线程说出了 hello，于是屏幕上会出现“hellohello world”。

\subsection{openmp: for 循环的并行}
用 openmp 中的 openmp for，可以将for循环并行化，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<omp.h>

int main(){

    #pragma omp parallel{//开启新的进程
        #pragma omp for //并行下面的for循环
        for(int i=1;i<=10;i++){
            cout<<"i="<<i<<"\t"<<"hello world!"<<endl;
        }
    }
    return 0;
}  
\end{lstlisting}
这个程序中的for循环被并行处理。
如果你的电脑有 8 个线程，for 循环有10次循环，那么每个线程自动分配得到其中的一部分任务（有的线程做1次循环，有的做2次循环）。第 9 行和第 11 行可以合并为一行：\#pragma omp parallel for。

如果不写第8行，那么每个线程都会执行一次9-11行，8个线程就会有80个“hello world”出现在屏幕上。

第7行和第8行可以合并为一行：\#pragma omp parallel for，所以上面的程序可以写作
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<omp.h>

int main(){

    #pragma omp parallel for //开启新的进程，并行下面的for循环
        for(int i=1;i<=10;i++){
            cout<<"i="<<i<<"\t"<<"hello world!"<<endl;
        }
    return 0;
}  
\end{lstlisting}
效果是一样的。

\subsection{线程的同步}
如果有两个线程同时读取1个公共变量，是没有问题的；但是如果有两个线程同时改变1个公共变量，就会起冲突。
有好几种办法避免这种冲突，\#pragma omp atomic 是其中的一种。
atomic 声明只对下方一行有效，这一行必须是如下形式之一：
\begin{itemize}
	\item x binop = expr
	\item x++/++x
	\item x- -/- -x
\end{itemize}
其中 x 是一个带值的变量（而不是数组），binop 可以是加减乘除运算符。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<omp.h>

int main(){

        double sum=0;
        #pragma omp parallel
        {
                #pragma omp for
                for(int i=1;i<=100;i++){
                        #pragma omp atomic  //保护下一行，任意时刻下一行只能被一个进程执行
                        sum += i;           //如果注掉#pragma omp atomic，1到100可能有些数漏加
                }
        }
        cout<<" Total sum from 1 to 100 = "<<sum<<endl;
        return 0;
}
\end{lstlisting}

如果要保护一整段代码，任一时刻只允许一个线程执行这段代码，则可以用 critical 命令，如下
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<omp.h>

int main(){

	double sum=0;
	#pragma omp parallel for 
	for(int i=1;i<=100;i++){
		#pragma omp critical
		{
			sum += i;
			if(i%10==0)cout<<"hello world"<<endl;//如果 i=10,20,...,100,则输出一个 hello world
		}
	}
	cout<<" Total sum from 1 to 100 = "<<sum<<endl;
	return 0;
}   
\end{lstlisting}
运行以后，你会看到，10个 hello world 整齐地出现，这是因为，任一时刻只允许一个线程执行这段代码，所以每个线程输出一个完整的 hello world 以后，下一个线程才会输出这句话。

atomic 设置付出的性能代价更小，所以，如果只有一行代码需要保护，就用 atomic 就行，不需要用 critical。

不同线程对任务的分配，可以有不同的策略，包括 static, dynamic 等等模式，可以查阅资料。
其中，static 是对运行时间影响最小的，也是系统默认的，所以在讲义里我们不予关注这个方面。

\subsection{openmp: reduction}
上面的高斯求和代码实现了并行，但实际上比串行程序还要慢，这是因为任意时刻，最多只能有一个进程执行
\begin{lstlisting}
sum += i;
\end{lstlisting}
所以所有进程都得排队，总的来说，程序还是在一个数一个数地，从１加到100。
为了改善这个情况，我们可以使用reduction。
\begin{lstlisting}
reduction(op:list)
\end{lstlisting}
op 表示运算符，包括$\{+,-,*\}$，list则表示一系列变量，用逗号隔开
\begin{itemize}
\item 每个进程会将list中的变量都copy一份局域变量，并赋初值(如果op为$+/-$,初值就是0，如果op是$*$，初值就是1)
\item 进程中的运算会更新这些局域变量
\item 这些局域变量按照op汇聚成一个最终结果，更新list中的变量。
\end{itemize}

下面是openmp编写的中值法积分代码
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<ctime>
#include<omp.h>
/*
 * integ_midpoint: 中值法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_midpoint(double a, double b, long long n, double (*func)(double)){

        double y=0;
        double h=(b-a)/n;
        #pragma omp parallel for reduction(+:y) // reduction
                for(long long i=0;i<n;i++){
                        y += func( a+ (i+0.5)*h);
                }
        y *= h;
        return y;
}

double f(double x){
        return x*x;
}

int main(){

        double a=0,b=1;
        long long n = 1E９;//长长整型，范围是[-2^{63}, 2^{63}-1]

        time_t t_start = time(NULL);
        cout<<"integ = "<<integ_midpoint(a,b,n,f)<<endl;
        time_t t_end = time(NULL);
        cout<<"Elapsed time: "<<(t_end-t_start)<<" s.\n";

        return 0;
}
\end{lstlisting}
在上面的代码中，integ\_midpoint()函数中实现了并行。
如果电脑有6个线程，则会赋值6个y的局域变量，并赋初值0，然后分给各个线程，每个线程负责自己的累加，都加到自己的y上，最后汇总，更新真正的y。

在主函数main中，我们用到了计时函数time，下面做简单讲解。
\paragraph{c++: 程序计时}
在ctime库中，time()函数与clock()函数都可以用来记录程序运行时间。
\vspace{1cm}
包括ctime库：
\begin{lstlisting}
#include<ctime>
\end{lstlisting}
使用time\_t类型，用time(NULL)得到当前时间
\begin{lstlisting}
time_t a = time(NULL);//从1970年1月1日0点整，到这一行执行的时刻，经过的总秒数
.....
time_t b = time(NULL);// 上面的程序段运行时间为　(b-a)s
\end{lstlisting}

使用clock\_t类型，用clock()得到{\bf cpu时间}。
\begin{lstlisting}
clock_t a = clock();//当前cpu时间，基本单位是　1/CLOCKS_PER_SEC　秒
.....
clock_t b = clock();// 上面的程序段运行的cpu时间为 (b-a)/CLOCKS_PER_SEC　秒
\end{lstlisting}
CLOCKS\_PER\_SEC在不同的系统上是不同的常数，在ubuntu中是$10^6$。

如果程序是单线程，这样用clock()记录的时间，差不多等于物理时间。
但如果程序是并行程序，如此得到的是所有线程的{\bf 总cpu时间}，可能并行以后程序变快了，但总cpu时间可能变多了。
一般来说，完成同一个任务，平均每线程的cpu时间只会随着线程数的增加而减小。
消耗的物理时间不会随着线程数的增大而无穷减小，而是很快饱和，在我有限的经验中，可能十来个线程就饱和了。

\section{OpenMPI：分布式内存并行}
因为时间关系，这部分内容不会在课程中讲授，只是放在这里供参考。
所以有些代码涉及到讲义后文中的内容，比如梯形法积分，当然这个也很简单。
\subsection{安装mpi, hello world!}
在ubuntu下非常简单
\begin{lstlisting}
sudo apt-get install mpich
\end{lstlisting}
下面是MPI并行的hello world程序
\begin{lstlisting}
#include<iostream>
using namespace std;
#include"mpi.h"
#include<cmath>

int main(int argc, char ** argv){

    int myid, numprocs;
    int namelen;
    char processor_name [ MPI_MAX_PROCESSOR_NAME ];

    MPI_Init( &argc, &argv);//每个进程都会复制上面的所有变量。
    MPI_Comm_rank( MPI_COMM_WORLD, & myid);//进程编号
    MPI_Comm_size( MPI_COMM_WORLD, & numprocs);//总进程数
    MPI_Get_processor_name( processor_name, & namelen);//主机名字

    fprintf( stderr, "Hello World! Process %d of %d on %s \n", myid, numprocs, processor_name);
    MPI_Finalize();
    return 0;
}
\end{lstlisting}
需要包含头文件"mpi.h"。其中用到了MPI中的几个函数和变量，
\begin{itemize}
\item MPI\_Init( int * argc\_p, char *** argv\_p ) 是初始化函数，做好 mpi 相关的初始化。传入的 argc\_p, argv\_p 是指向 main 函数的命令行参数 argc, argv 的指针。如果没有命令行参数，也可以输入 NULL。

\item MPI\_COMM\_WORLD：通讯子，“一组可以互发消息的进程集合”

\item int MPI\_Comm\_rank( \_\_in MPI\_Comm comm, \_\_out int * rank ) ：第一个参数是输入参数——通信子，第二个参数是输出参数——进程号。

\item int MPI\_Comm\_size( \_\_in MPI\_Comm comm, \_\_out int * size )：第一个参数是输入参数——通信子，第二个参数是输出参数：通信子中总进程数。

\item MPI\_MAX\_PROCESSOR\_NAME：似乎是一个整数，值为128.

\item MPI\_Get\_processor\_name( ..)：得到主机的名字

\paragraph{编译，运行！}
编译MPI并行的cpp程序：
\begin{lstlisting}
mpicxx main.cpp
\end{lstlisting}
在linux系统下，会自动生成可执行文件a.out.
运行可执行文件：
\begin{lstlisting}
mpirun -np 2 ./a.out
\end{lstlisting}
np表示number of processes，即并行的进程数，所以上面的命令用2个进程执行并行任务。
运行的结果为
\begin{lstlisting}
Hello World! Process 0 of 2 on Swagger 
Hello World! Process 1 of 2 on Swagger
\end{lstlisting}
即每个进程都输出自己的进程编号、总进程数、以及电脑名字。

如果在MPI\_Finalize后面再写几行内容，这几行也会被各个进程都执行一遍。
所以，似乎MPI的程序没法进行串行->并行->串行的计算，只能从头到尾都是并行。

\subsection{点对点通信：MPI\_Send和MPI\_Recv}
两个特定进程之间进行通信，就叫做点对点通信，可以用MPI\_Send和MPI\_Recv实现。
下面的代码将所有进程编号都加起来。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include"mpi.h"
#include<cmath>

int main(int argc, char ** argv){

    int myid, numprocs;
    int namelen;
    char processor_name [ MPI_MAX_PROCESSOR_NAME ];
    int ha;
    int i;
    int total = 0;

    MPI_Init( &argc, &argv);
    MPI_Comm_rank( MPI_COMM_WORLD, & myid);
    MPI_Comm_size( MPI_COMM_WORLD, & numprocs);


    if( myid !=0 ){
        ha = myid;//进程编号
        MPI_Send( &ha, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);//将1整型ha的值向进程0发送
    }
    else{
        for( i = 1; i < numprocs; i++){
            MPI_Recv( &ha, 1, MPI_INT, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            //从进程i接收1个整型ha的值
            total += ha;
        }
        cout<< " total = "<<total<<endl;
    }

    MPI_Finalize();

    return 0;
}
\end{lstlisting}
\item int MPI\_Send (void *buf, int count, MPI\_Datatype datatype,int dest, int tag,MPI\_Comm comm)：点对点发送消息。
    \begin{itemize}
    \item *buf: 发送缓冲区，指向发送的变量地址
    \item count: 发送的数据个数
    \item datatype：发送的数据类型
    \item dest: 目标进程序号，取值为 0 到 np-1 的整数
    \item tag: 消息标签，取值为 0  到 MPI\_TAG\_UB
    \item comm：通信子，即线程集合
    \end{itemize}

\item int MPI\_Recv (void *buf, int count, MPI\_Datatype datatype,int source, int tag, MPI\_Comm comm,MPI\_Status *status)：点对点接收消息。与 MPI\_Send 相似，不同的是 MPI\_Status * status 为接收状态。MPI\_Status 似乎是 MPI 自己定义的一个类。
\end{itemize}

\subsection{集体通信：MPI\_Reduce}
进程之间的通信方式为1对多、多对1、或者多对多，就称作是集体通信。
下面的代码用MPI\_Reduce将所有进程编号加起来。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include"mpi.h"
#include<cmath>

int main(int argc, char ** argv){

    int myid, numprocs;
    int namelen;
    char processor_name [ MPI_MAX_PROCESSOR_NAME ];
    int ha;
    int i;
    int total = 0;

    MPI_Init( &argc, &argv);//初始化
    MPI_Comm_rank( MPI_COMM_WORLD, & myid);//进程编号
    MPI_Comm_size( MPI_COMM_WORLD, & numprocs);//总进程数

    MPI_Reduce( & myid, & total, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
    //所有进程都将自己的1个整型数据——myid值——发送给进程0的total数据，汇聚方式为累加。

    if( myid ==0 ){//进程0输出total的结果
        cout<<" total = "<<total<<endl;
    }

    MPI_Finalize();

    return 0;
}
\end{lstlisting}
MPI\_Reduce的函数头如下
\begin{lstlisting}
 int MPI_Reduce (void *sendbuf, void *recvbuf, int count,MPI_Datatype datatype, MPI_Op op, int root,MPI_Comm comm)
\end{lstlisting}
各参数意义是自明的。

下面我们举几个例子，展示MPI的简单应用。
\subsection{MPI：梯形法数值积分}
梯形法数值积分的公式为
\begin{eqnarray}
I_f &\approx& \sum^{n-1}_{i=0} h (f(x_i) + f(x_{i+1}))/2 \nonumber\\
    &=& h \{ \frac{f(x_0)}{2} + f(x_1) + \cdots + f(x_{n-1}) + \frac{f(x_n)}{2} \}.
\end{eqnarray}
下面的代码将积分区域分割，将积分工作分派给不同的进程，然后用MPI\_Send和MPI\_Recv通信，收集起各个进程的结果，汇总得到整个积分值。
\begin{lstlisting}
#include<iostream>
using namespace std;

#include"mpi.h"
#include<cmath>

double f(double x){
    return pow(x,3);
}

double Trap(double left_endpt, double right_endpt, double trap_count, double base_len){

    double estimate = ( f(left_endpt) + f(right_endpt) )/2, x = left_endpt;

    for(int i=1;i<=trap_count-1;i++){
        x += base_len;
        estimate += f(x);
    }
    estimate = estimate * base_len;
    return estimate;
}

int main(){

    int my_rank = 0, comm_sz = 0, n = 1E7, local_n = 0;
    double a = 0, b = 3, h = 0, local_a = 0, local_b = 0;
    double local_int = 0, total_int = 0;
    int source;
    clock_t t_start = clock();

    MPI_Init( NULL, NULL );
    MPI_Comm_rank( MPI_COMM_WORLD, & my_rank );
    MPI_Comm_size( MPI_COMM_WORLD, & comm_sz );

    double start_time = MPI_Wtime();

    h = (b-a)/n; 
    local_n = n / comm_sz;
    if( my_rank == comm_sz - 1 ){
        local_n = n - local_n * (comm_sz-1);
    }// when n is not indivisible by comm_sz, the last chunk is larger than others.

    local_a = a + my_rank * ( n / comm_sz ) * h;
    local_b = local_a + local_n * h;
    local_int = Trap( local_a, local_b, local_n, h );
    //printf("%d of %d processes, local_int = %.15e\n", my_rank, comm_sz, local_int);

    if( my_rank != 0 ){
        MPI_Send( & local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD );
    }
    else{
        total_int = local_int;
        for( source = 1; source < comm_sz; source ++ ){
            MPI_Recv( & local_int, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE );
            total_int += local_int;
        }

        printf("With n = %d trapezoids, our estimate \n", n);
        printf("of the integral from %f to %f = %.15e\n", a, b, total_int );

        double end_time = MPI_Wtime();
        printf("That took %lf seconds.\n", end_time - start_time);    
    }

    MPI_Finalize();

    return 0;
}
\end{lstlisting}
其中用到了MPI\_Wtime()函数，这个函数没有参数，它选一个过去的时刻（具有任意性），然后计算调用时距离那个选定时刻的物理时间差，所以它的返还值意义不大，但两次调用得到的时间差可以用来给程序计时。
MPI\_WTICK决定了返还值的精度，一般返还值精度很高。

图\ref{fig:MPI-trapzoid}是进程数不同时，程序的耗时。
\begin{figure} [h]
\centering
\includegraphics[width=\textwidth]{MPI_trapzoid}
\caption{进程数不同时，MPI梯形法积分的耗时。}
\label{fig:MPI-trapzoid}
\end{figure}

\subsection{矩阵相乘}
下面的代码用MPI\_Recv和MPI\_Send实现矩阵分块相乘。
为了计算矩阵C=A*B，进程0将A分为排成一列的若干分块矩阵，然后将这些分块矩阵与矩阵B分发给不同进程，这些进程计算得到C的分块矩阵，发送给进程0，进程0汇总，完成任务。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<mpi.h>

int main(int argc, char * argv[] ){

    double start, stop;
    int *a, *b, *c, *buffer, *ans;
    int size = 1000;
    int rank, numprocs, line;

    MPI_Init(NULL,NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    line = size/numprocs;

    b = new int [ size * size ];
    ans = new int [ size * line ];

    start = MPI_Wtime();

    if( rank ==0 ){

        a = new int [ size * size ];
        c = new int [ size * size ];

        for(int i=0;i<size; i++)
        for(int j=0;j<size; j++){
            a[ i*size + j ] = i*j;
            b[ i*size + j ] = i + j;
        }

        for(int i=1;i<numprocs;i++){// send b
            MPI_Send( b, size*size, MPI_INT, i, 0, MPI_COMM_WORLD );
        }

        for(int i=1;i<numprocs;i++){// send part of a
            MPI_Send( a + (i-1)*line*size, size*line, MPI_INT, i, 1, MPI_COMM_WORLD);
        }

        for(int i = (numprocs-1)*line;i<size;i++){// calculate block 1
            for(int j=0;j<size;j++){
                int temp = 0;
                for(int k=0;k<size;k++)
                    temp += a[i*size+k]*b[k*size+j];
                c[i*size+j] = temp;
            }
        }

        for(int k=1;k<numprocs;k++){// recieve ans
            MPI_Recv( ans, line*size, MPI_INT, k, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            for(int i=0;i<line;i++){
                for(int j=0;j<size;j++){
                    c[ ((k-1)*line + i)*size + j] = ans[i*size+j];
                }
            }
        }

        FILE *fp = fopen("c.txt","w");
        for(int i=0;i<size;i++){
            for(int j=0;j<size;j++)
                fprintf(fp,"%d\t",c[i*size+j]);
            fputc('\n',fp);
        }
        fclose(fp);
        
        stop = MPI_Wtime();

        printf("rank:%d time:%lfs\n",rank,stop-start);

        delete [] a,c;

    }
    else{
        buffer = new int [ size * line ];

        MPI_Recv(b, size*size, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(buffer, size*line, MPI_INT, 0, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        for(int i=0;i<line;i++)
        for(int j=0;j<size;j++){
            int temp=0;
            for(int k=0;k<size;k++)
                temp += buffer[i*size+k]*b[k*size+j];
            ans[i*size+j] = temp;
        }
        MPI_Send(ans, line*size, MPI_INT, 0, 3, MPI_COMM_WORLD);

        delete [] buffer;
        delete [] ans;
    }

    delete [] b;

    MPI_Finalize();

    return 0;
}
\end{lstlisting}

也可以用集体通信，完成同样的任务，下面的代码用MPI\_Scatter和MPI\_Gather，完成矩阵相乘C=A*B。
MPI\_Scatter从一个进程向多个进程发送消息，它的函数头为（参数意义是自明的），
\begin{lstlisting}
int MPI_Scatter(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
               void *recvbuf, int recvcount, MPI_Datatype recvtype, int root,
               MPI_Comm comm)
\end{lstlisting}


MPI\_Gather从多个进程收集数据，它的函数头为
\begin{lstlisting}
int MPI_Gather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
               void *recvbuf, int recvcount, MPI_Datatype recvtype,
               int root, MPI_Comm comm)
\end{lstlisting}

代码如下
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<mpi.h>

int main(){

    int my_rank;
    int num_procs;
    int size = 1000;
    double start, finish;

    MPI_Init(NULL,NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    
    int line = size / num_procs;
    cout<<" line = "<<line<<endl;
    int * local_a = new int [ line * size ];
    int * b = new int [ size * size ];
    int * ans = new int [ line * size ];
    int * a = new int [ size * size ];
    int * c = new int [ size * size ];

    if( my_rank == 0 ){

        start = MPI_Wtime();

        for(int i=0;i<size;i++){
            for(int j=0;j<size;j++){
                a[ i*size + j ] = i*j;
                b[ i*size + j ] = i + j;
            }
        }

        MPI_Scatter(a, line * size, MPI_INT, local_a, line * size, MPI_INT, 0, MPI_COMM_WORLD );

        MPI_Bcast(b, size*size, MPI_INT, 0, MPI_COMM_WORLD);

        for(int i= 0; i< line;i++){
            for(int j=0;j<size;j++){
                int temp = 0;
                for(int k=0;k<size;k++)
                    temp += a[i*size+k] * b[k*size + j];
                ans[i*size + j ] = temp;
            }
        }

        MPI_Gather( ans, line * size, MPI_INT, c, line * size, MPI_INT, 0, MPI_COMM_WORLD );

        for(int i= num_procs *line; i< size;i++){
            for(int j=0;j<size;j++){
                int temp = 0;
                for(int k=0;k<size;k++)
                    temp += a[i*size+k] * b[k*size + j];
                c[i*size + j ] = temp;
            }
        }

        FILE *fp = fopen("c2.txt","w");
                for(int i=0;i<size;i++){
                        for(int j=0;j<size;j++)
                                fprintf(fp,"%d\t",c[i*size+j]);
                        fputc('\n',fp);
                }
                fclose(fp);

        finish = MPI_Wtime();
        printf(" time: %lf s \n", finish - start );
    }
    else{
        int * buffer = new int [ size * line ];
        MPI_Scatter(a, line * size, MPI_INT, buffer, line * size, MPI_INT, 0, MPI_COMM_WORLD );
        MPI_Bcast( b, size * size, MPI_INT, 0, MPI_COMM_WORLD );
        /*
        cout<<" b:"<<endl;
        for(int i=0;i<size;i++){
            for(int j=0;j<size;j++){
                cout<<b[i*size + j]<<",";
            }
            cout<<endl;
        }
        */

        for(int i=0;i<line;i++){
            for(int j=0;j<size;j++){
                int temp = 0;
                for(int k=0;k<size;k++)
                    temp += buffer[i*size+k] * b[k*size + j];
                //cout<<"i = "<<i<<"\t j= "<<j<<"\t temp = "<<temp<<endl;
                ans[i*size + j] = temp;
            }
        }

        /*
        cout<<" ans:"<<endl;
        for(int i=0;i<line;i++){
            for(int j=0;j<size;j++){
                cout<<ans[i*size + j]<<",";
            }
            cout<<endl;
        }
        */

        MPI_Gather(ans, line*size, MPI_INT, c, line*size, MPI_INT, 0, MPI_COMM_WORLD );
        delete [] buffer;
    }

    delete [] a, local_a, b, ans, c;

    MPI_Finalize();
    return 0;
}
\end{lstlisting}

下面是这两种方式做矩阵相乘的消耗时间，
\begin{figure} [h]
\centering
\includegraphics[width=\textwidth]{MPI_matrix_multiplication}
\caption{进程数不同时，MPI矩阵相乘的耗时。}
\label{fig:MPI-matrix-multiplication}
\end{figure}

MPI还有很多库函数，我给出官网中函数索引的链接\footnote{http://www.mpich.org/static/docs/v3.2/www3/}。

\section{Linux系统：Ubuntu}
在这一节中，我放进一些关于linux系统ubuntu18.04的笔记。
这些内容不会在课堂上讲，仅仅供参考。

\subsection{rar：处理winrar文件}
Ubuntu下可以安装rar，
\begin{lstlisting}
sudo apt-get install rar
\end{lstlisting}
如果要新建a.rar文件，将a.out压缩进去
\begin{lstlisting}
rar a a.rar a.out
\end{lstlisting}
如果要解压缩a.rar文件
\begin{lstlisting}
rar x a.out
\end{lstlisting}

\subsection{vim中文编码}
vim是常用的工具，但如果文件中有gbk等编码，而vim默认编码是utf8，中文部分可能就成为乱码。
可以做如下设置，解决这个问题：用root权限打开/etc/vim/vimrc文件，在文件末尾加上
\begin{lstlisting}
   set fileencodings=utf-8,ucs-bom,gb18030,gbk,gb2312,cp936
   set termencoding=utf-8
   set encoding=utf-8
\end{lstlisting}
设定了文件的解码选项，其中包括了gbk相关的编码，所以就可以正常显示中文了。

\subsection{命令行发送邮件}
安装 sendemail
\begin{lstlisting}
	sudo apt-get install sendemail
\end{lstlisting}
取得你邮箱地址的 stmp 授权码，比如 126 邮箱，登录邮箱以后，找到设置->POP3/SMTPIMAP，然后按提示开通 STMP 服务，设置 STMP 授权码。
有了这个授权码，你才有权限在 linux 命令行下用这个邮箱发送邮件。
然后在 Ubuntu 下，创建一个文件，比如
\begin{lstlisting}
	touch /home/luyi/Desktop/hello.sh
\end{lstlisting}
即创建一个脚本，打开脚本
\begin{lstlisting}
	vi hello.sh
\end{lstlisting}
在其中输入
\begin{lstlisting}
sendEmail -xu luyi19881007@126.com -xp smtp授权码 -t 目标邮箱 -u "Greetings" -m "Have a good day!" -s smtp.126.com -f luyi19881007@126.com
\end{lstlisting}
其中 -xu 后面是你的邮箱地址，-xp 后面是你的邮箱 stmp 授权码，-t 后面是目标邮箱地址，-u 后面双引号内是主题，-m 后面双引号内是邮件内容，-s 后面是 126 邮箱的 stmp 域名，如果是 qq 邮箱，需要查阅，-f 后面是发送邮箱地址。
在上本课程期间，大家将在节假日收到课程邮箱的祝福。

写了这个脚本以后，每次发送邮件，只需执行该脚本
\begin{lstlisting}
	sh /home/luyi/Desktop/hello.sh
\end{lstlisting}
系统就会提示，邮件是否发送成功。

\subsection{ubuntu系统定时执行任务}
我们要让系统自动执行上面的脚本，即可自动发送邮件出去。
那么怎么让系统定时执行某项任务呢？
ubuntu 系统自带 cron 命令，可以设置定时执行的任务。
我们需要设置 cron 的配置文件，让它在相应时间做相应的事。
打开 /etc/crontab 文件，
\begin{lstlisting}
	sudo vim /etc/crontab
\end{lstlisting}
在文件的最下方加上一行
\begin{lstlisting}
0 8   * * *   root    sh /home/luyi/Desktop/hello.sh
\end{lstlisting}
这一行的各个数字的意义分别是：\\
m:			minute\\
h:			hour\\
dom: 		day of month\\
mon:		month of year\\
dow:		day of week\\
user:		用户\\
command:	执行的任务\\
所以我加的那一行的意思是：“每一天早上 8 点整执行 hello.sh 这个脚本”。

\subsection{ubuntu：开机自动启动脚本}
我们最后还需要做一步，就大功告成了。
cron 这个工具并不会开机自动启动，所以我们需要手动设置。
打开 /etc/rc.local 文件（如果没有，可以自己手动添加）,并加入一行
\begin{lstlisting}
service cron start

exit 0
\end{lstlisting}
其中，service 是一个 命令，可以启动、关闭、重启一个任务，所以这一行的意思是：“启动 cron”。
exit 0 的意思是正常退出。
保存并退出，就大功告成了。

这样，只要电脑开着，就会自动按指定要求发送邮件。

\subsection{linux：脚本编写}
在linux系统下，有时一段命令会重复使用，每次都重敲一遍很浪费时间。。。
有时需要让电脑依次运行一系列任务，每个任务都需要一段时间，人类等着发号施令，也很无聊。。。
所以，可以把命令写到一个文本文件里，让电脑自动依次执行其中的命令，这样的文本文件就叫做脚本。
如下是一个hello world脚本
\begin{lstlisting}
#!/bin/bash
echo "Hello World!"
\end{lstlisting}
第一行是注释，注释行以\#开头。
注释一般是不执行的，只起说明的作用，但如果是第一行，并且以“\#!”开头，就是唯一的例外：上面的例子中第一行告诉系统，用/bin/bash执行这个脚本。
也就是说，把上面的脚本存为hello.sh，然后赋予脚本可执行权，
\begin{lstlisting}
chmod u+x hello.sh
\end{lstlisting}
然后
\begin{lstlisting}
./hello.sh
\end{lstlisting}
会得到
\begin{lstlisting}
Hello World!
\end{lstlisting}

\paragraph{if语句}
每一个命令行执行后会有一个返还值，可以用\$?查看，比如：
\begin{lstlisting}
echo $?
\end{lstlisting}
如果上一次命令正常执行了，会得到
\begin{lstlisting}
0
\end{lstlisting}
否则，可能得到1-254间的任何一个整数。

利用命令返还的值，可以写if语句，根据情况输出不同的结论
\begin{lstlisting}
#!/bin/bash
# testing nested ifs
#
#
testuser=NoSuchUser
#
if grep $testuser /etc/passwd
then
        echo "The user $testuser exists on this system."
#
elif ls -d /home/$testuser
then
        echo "The user $testuser does not exist on this system."
        echo "However, $testuser has a directory."
#
else
        echo "The user $testuser does not exist on this system."
        echo "And, $testuser does not have a directory."
fi
\end{lstlisting}
上面的代码判断 /etc/passwd 中是否有某个用户的文件，并进一步判断 /home/\$testuser 文件夹是否存在，根据结果输出相应结论。

\paragraph{case语句}
case语句的语法为
\begin{shaded*}
\noindent case variable in \\
pattern1 | pattern2) commands1;; \\
pattern3 ) commands2;; \\
*) commands3;; \\
esac
\end{shaded*}
例子：
\begin{lstlisting}
#!/bin/bash

case $USER in
luyi | haha | hehe )
        echo "Welcome, $User"
        echo "Please enjoy your visit";;
luhan )
        echo "Please don't forget to logout, thx.";;
*)
        echo "Sorry, you are not allowed here.";;
esac
\end{lstlisting}

\paragraph{for语句}
for语句的语法为
\begin{shaded*}
\noindent for var in list \\
do \\
\indent    commands \\
done 
\end{shaded*}
例如
\begin{lstlisting}
for test in Alabama Alaska Arizona Arkansas California Colorado
do
    echo The next is $test
done
\end{lstlisting}
c语言风格的for循环
\begin{lstlisting}
#!/bin/bash
# c style for
echo "example for C style for:"
for((a=1;a<10;a++))
do
        echo "a = $a"
done

#multiple variables
echo
echo "example for multiple variables in C style for:"
for((a=1,b=10; a<=10; a++,b--))
do
        echo "$a - $b"
done

#nested 
echo
echo "example for nested C style for:"
for((a=1;a<=3;a++))
do
        echo "Starting loop $a:"
        for((b=1;b<=3;b++))
        do
                echo " Inside loop: $b"
        done
done
\end{lstlisting}
上面的代码中，展示了c语言风格的for循环、多变量的for循环，以及嵌套的for循环。

\paragraph{换行符}
如果一行命令太长，需要用几行来表示，那么可以在换行出用$\backslash$，即表示命令没有完。
但是需要注意的是，如果一个命令行分为两部分，中间被一行注释隔开，是会被当成两行命令的。

\subsection{远程登录：ssh}
如果ssh服务端的端口为1105，ip为202.120.52.21，用户名为Daddy，则如下登录服务器
\begin{lstlisting}
ssh -p1105 Daddy@202.120.52.21
\end{lstlisting}
如果端口、用户名、ip都正确，服务器会要求输入用户名的密码，若输入正确，就会成功登录。

如果要远程传送文件，
\begin{lstlisting}
ssh -P1105 ha.txt Daddy@202.120.52.21:~/
\end{lstlisting}

\paragraph{自动输入密码：sshpass}
有时候需要写脚本，自动向服务器上传文件，需要自动输入密码，就可以安装sshpass
\begin{lstlisting}
sudo apt-get install sshpass
\end{lstlisting}
然后如下登录/传送文件
\begin{lstlisting}
sshpass -p passwd ssh -p1105 Daddy@202.120.52.21
sshpass -p passwd scp -P1105 ha.txt Daddy@202.120.52.21:~/
\end{lstlisting}
也就是说，在原来的ssh/scp命令前面加上sshpass -p passwd即可，passwd是相应的密码。

\subsection{linux:云服务器}
这部分内容是我的师兄雷杨教我的，哈哈哈。
\paragraph{申请amazon云服务器}
\begin{itemize}
\item 打开网址　https://aws.amazon.com/cn/free
\item 按照提示，申请一个免费云服务器，选择 Ubuntu18.04，并启动服务器。中间会下载一个private key文件，保存到一个安全的位置，我保存在了/etc/network/下面，这样我就不会误删它。
\end{itemize}

\paragraph{云服务器作为代理：firefox设置}
\begin{itemize}
\item Preferences -> Network Settings -> Manuel proxy settings -> SOCKS v5
\item SOCKS Host: 127.0.0.1, port: 1080, 
\item 地址栏输入：about:config，选择 I accept the risk!，在搜索栏输入 dns，找到network.proxy.socks\_remote\_dns，双击变成true即可
\end{itemize}

\paragraph{云服务器作为代理：登录！}

打开terminal, 输入
\begin{lstlisting}
ssh -D 1080 -i your_private_key.pem ubuntu@your_server_address
\end{lstlisting}
然后用firefox浏览，即可。
